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Abstract. We present new results from accurate and fully general-relativistic simulations of 
the coalescence of unmagnetized binary neutron stars with various mass ratios. The evolution 
of the stars is followed through the inspiral phase, the merger and prompt collapse to a 
black hole, up until the appeai'ance of a thick accretion disk, which is studied as it enters 
and remains in a regime of quasi-steady accretion. Although a simple ideal-fluid equation 
of state with F = 2 is used, this work presents a systematic study within a fully general 
relativistic framework of the properties of the resulting black-hole-torus system produced 
by the merger of unequal-mass binaries. More specifically, we show that: (1) The mass of 
the torus increases considerably with the mass asymmetry and equal-mass binaries do not 
produce significant tori if they have a total baryonic mass A/tot ^ 3.7 Mq; (2) Tori with 
masses Mtor ~ 0.2 Mq are measured for binaries with A/tot ~ 3.4 Mq and mass ratios 
q ~ 0.75 — 0.85; (3) The mass of the torus can be estimated by the simple expression 
Mtor(<7, Mtot) = [ci(l - <j) + C2] (Mniax — Mtot), involving the maximum mass for the 
binaries and coefficients constrained from the simulations, and suggesting that the tori can 
have masses as large as Aftor ~ 0.35 Mq for Mtot ~ 2.8 Mq and q ~ 0.75 — 0.85; 
(4) Using a novel technique to analyze the evolution of the tori we find no evidence for the 
onset of non-axisymmetiic instabihties and that very Httle, if any, of their mass is unbound; (5) 
Finally, for all the binaries considered we compute the complete gravitational waveforms and 
the recoils imparted to the black holes, discussing the prospects of detection of these sources 
for a number of present and future detectors. 



PACS numbers: 04.30.Db, 04.40.Dg, 04.70.Bw, 95.30.Lz, 97.60.Jd 

1. Introduction 

The numerical investigation of the coalescence and merger of binary neutron stars (NSs) 
within the framework of general relativity is receiving increasing attention in recent years 
(see e.g. |[Tl|2][3l|4l|5]|6ll2l[ll and references therein). Drastic improvements in the simulation 
front regarding mathematics {e.g. formulation of the equations), physics (e.g. incorporation 
of equations of state (EOSs) from nuclear physics), and numerical methods (e.g. use of 
high-resolution methods and adaptive mesh refinement), along with increased computational 
resources, have allowed to extend the scope of the early simulations (e.g. |9l). Larger initial 
separations have recently started being considered and some of the existing simulations have 
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expanded the range spanned by the models well beyond black-hole (BH) formation [B] [8] |71 . 
This is allowing the computation of the entire gravitational waveform from the early inspiral 
to the decaying tail of the late ringing of the formed BH. The construction of such waveform 
templates is still one of the driving motivations to perform binary NS simulations, as such 
events are among the most promising sources of detectable gravitational radiation for laser 
interferometric detectors. The current estimate for the detection rate relative to the first- 
generation interferometric detectors is ~ 1 event per ^ 40 — 300 years, increasing to an 
encouraging ^ 10 — 100 events per year for the advanced detectors lITOl . The second major 
incentive behind this type of simulations is establishing whether the end-product of the merger 
can act as the underlying mechanism operating at the central engine of short-hard gamma-ray 
bursts (SGRBs) lfm[T2ll . The consensus emerging from the existing simulations suggests the 
formation, depending on the suitability of the initial parameters of the simulated model, of 
a BH of stellar mass surrounded by a hot disk. Driven by neutrino processes and magnetic 
fields, such a compact system may be capable of launching a relativistic fireball with an energy 
of ^ 10"** erg on a timescale of 0.1 — Is ifTSl . 

This paper is dedicated in particular to investigating the late-time dynamics of the torus 
formed after the merger of unequal-mass NS binaries. As we describe below, all but one model 
of our initial sample have mass ratio different from one. Previous simulations have shown that 
the key parameter controlling the amount of mass left in the disk for a given initial mass in the 
system and EOS is the NS mass ratio IIT4l [n. Broadly speaking the general trend is simple: 
the larger the departure from equal-mass ratio, the more important tidal effects become in the 
less massive star, resulting in its tidal disruption. Because this takes place when the separation 
is still comparatively large, the angular momentum of the matter is still large and it results in 
larger-size and more massive disks. Early and low-resolution simulations with an ideal-gas 
EOS lfT4l have been shown to yield a disk mass of several percents of the total mass of the 
system for a mass ratio of ^ 0.85. Improved simulations by IJ] which adopted a hybrid 
EOS to mimic realistic, stiff nuclear EOS, indicate that the mass of the disk is ^ 0.01 Mq or 
slightly larger when the merger does not result in prompt collapse to a BH but in the formation 
of a hypermassive NS of large ellipticity (which later collapses to a BH following angular- 
momentum transport caused by emission of gravitational radiation). Similar disk masses, as 
large as ~ 0.02 Mq, are also reported in the latest simulations of |[8l, in which the initial 
orbital separation of the two stars is larger than in previous works. 

Observational data seems to indicate that the total gravitational masses of the known 
galactic NS binary systems are in a narrow range ~ 2.65 — 2.85 Mq ifTSl and there is also 
evidence indicating that the masses of the two NSs are nearly equal, with the baryonic mass 
ratio q = M1/M2 being between 1 and ^ 0.7 (Hereafter we will refer to q simply as the 
"mass ratio" but we report in table [T] also the ratio in the ADM masses (Z^dm ' 1 9adm 
do not coincide because of the nonlinear relation between baryonic and gravitational mass). 
However, there is no theoretical reason to assume that unequal-mass NS binaries are not 
produced in nature as often as the seemingly prevailing almost exactly equal-mass systems, 
particularly for twin giant progenitors IfTHfTTI . There are, indeed, recent computations and 
observations that contradict the predominance of almost exactly equal-mass systems ifTSlfTTl . 
On the one hand, binary population synthesis computations performed by 1 1 8!] show two peaks 
in the observability-weighted distribution of double NSs. One of these peaks is around q ^ \ 
and appears when both masses are close to 1.4 Mq. The second peak is around considerably 
smaller mass ratios and depends on the assumed maximum mass of a single NS (which is in 
turn dependent on the EOS considered): the higher this mass (the stiffer the EOS) the more 
significant the second peak ifTSl |T9l . However, the crucial parameter determining the shape 
of the distribution is the inclusion of hypercritical accretion onto the compact object during 
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the brief but critical "common-envelope" evolution phase of the close binary 119] . Similarly, 
recent computations by ifTTl also accounting for the effects of hypercritical accretion during 
the red-giant evolution of the less massive component of the binary lead to a pattern of NS 
binaries consisting of pulsars which are ^ 50% more massive than their companion NSs. 

An additional issue which motivates our study has to do with the investigation of 
the long-term stability and dynamics of the formed accretion disks. It is well known 
that thick accretion disks orbiting BHs may be subject to a number of instabilities, both 
axisymmetric, such as the so-called "runaway instability" ll20l . or non-axisymmetric, such 
as the "Papaloizou-Pringle instability" ||2T]| . The first one, in particular, if present, could 
destroy the torus on dynamical timescales, challenging the viability of the BH-torus 
model for SGRBs. Early time-dependent, axisymmetric, general-relativistic-hydrodynamical 
simulations of the runaway instability of non-self-gravitating tori around BHs were performed 
by II22I I23I . The distribution of specific angular momentum in the disk, £ = —u^/ut, 
with it0 and ut being the corresponding components of the 4-velocity u^, was the key 
parameter discriminating stable from unstable models in those simulations. It was found 
that ^—constant models were runaway unstable while power-law distributions I = iiTr" 
were stable for very small values of the angular momentum slope a (much smaller than 
the Keplerian value a = 0.5). Recent fully relativistic simulations by ll24ll . which for 
the first time have take into account the self-gravity of the disk, indicate that self-gravity 
does not favour the appearance of the instability, irrespective of the angular-momentum 
distribution. Under the effect of a perturbation, marginally stable models show the presence 
of axisymmetric oscillations for several dynamical timescales without the manifestation of 
the runaway instability, as 1251 had previously found for the case of non-self-gravitating tori. 
Indeed, the introduction of perturbations triggers quasi-periodic oscillations (QPOs) lasting 
tens of orbital periods, with amplitudes that are modified only slightly by the small loss of 
matter across the cusp 12511261 . The spectral distribution of the associated eigenfrequencies 
shows the presence of a fundamental p mode and of a series of overtones in a harmonic ratio 
2:3, which have been proposed to explain the QPOs observed in the X-ray luminosity of low- 
mass X-ray binaries (LMXBs) containing a BH candidate with the QPOs of small tori near 
the BH l27l l28l |29l . In addition, when sufficiently massive and compact, the oscillations of 
these tori are responsible for an intense emission of gravitational waves l25l [30l [311 l26ll . 

Overall, the ab-initio simulations reported here indicate that large-scale tori with masses 
Mtor ~ 0.2 Mq can be produced as the result of the inspiral and merger of binary NSs with 
unequal masses and that even larger masses can be predicted for binaries with smaller total 
masses. These tori are typically of large size, with quasi-Keplerian distribution of angular 
momenta, showing quasi-stationary evolutions and the absence of dynamical instabilities. As 
such, these results may provide additional information relevant to all the above issues for BH- 
torus systems formed in a fully consistent manner within the framework of general relativity. 
Furthermore, the gravitational-wave emission computed here reveals that the waveforms are 
sensitive to the mass ratio in the binary, both during the inspiral and after the merger, and 
could be used to extract important information on the structure and EOS of the progenitor 
stars. Such observations, however, will most likely have to rely on the advanced detectors, 
which will become operative in a few years. 

The paper is organised as follows; Section |2] describes the mathematical and numerical 
framework of our simulations. Section |3] discusses the dynamics of the coalescence and 
merger of our model sample. Next, in Section |4] we focus on the analysis of the tori 
formed after the merger and on their physical properties. The issue of the gravitational- 
wave emission from unequal-mass NS mergers is discussed in Section |5] and the main 
conclusions of our investigation are presented in Section|6] In addition Appendix A provides 
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quantitative measures of the accuracy of our numerical methods. We use a spaceUke signature 
(—,+,+,+) and a system of units in which c = G = Mq = 1 (or in other units whenever 
more convenient). Greek indices are taken to run from to 3, Latin indices from 1 to 3 and 
we adopt the standard convention for the summation over repeated indices. 

2. Mathematical and Numerical Setup 

All the details on the mathematical and numerical setup used for producing the results 
presented here are discussed in depth in 1321 [31. In what follows, we limit ourselves to a 
brief overview. 

2.1. Einstein and Hydrodynamics equations 

The evolution of the spacetime is obtained using the CCATIE code, a three-dimensional 
finite-differencing code providing a solution of a conformal traceless formulation of the 
Einstein equations with a "1 + log" slicing condition and a "Gamma-driver" shift condition 
(the interested reader is addressed to 1321 [33ll for a detailed discussion of the equations and 
gauges used). The general-relativistic hydrodynamics equations are instead solved using the 
Whisky code presented in l34l [3511331 . which adopts a flux-conservative formulation of the 
equations as presented in l36ll and high-resolution shock-capturing schemes. The Whisky 
code implements several reconstruction methods, such as Total- Variation-Diminishing (TVD) 
methods, Essentially-Non-Oscillatory (ENO) methods l37l . and the Piecewise Parabolic 
Method (PPM) l38l . Also, a variety of approximate Riemann solvers can be used, starting 
from the Harten-Lax-van Leer-Einfeldt (HLLE) solver l39l , over to the Roe solver l40ll . and 
the Marquina flux formula l4Tll (see l34ll35ll for a more detailed discussion). All the results 
reported hereafter have been computed using the Marquina flux formula 1421 and a PPM 
reconstruction. We stress again (as already done in l3]|2l) that the use of high-order methods 
and high-resolution is essential to be able to draw robust conclusions on the inspiral and 
merger. Lower-order methods in the reconstruction and low resolution may yield convergent 
and apparently reasonable results which however contain a large truncation error Specific 
examples of this type of problem are presented in Appendix A of Q and in Figure 4 of Q. 
Also, a measure of our overall accuracy is presented in Appendix A below and shows that by 
employing such methods we are able to conserve energy and angular momentum to ~ 1% 
over a timescale of ~ 140 ms. 

The system of hydrodynamics equations is closed by an EOS. As discussed in detail 
in 121, the choice of the EOS plays a fundamental role in the post-merger dynamics and 
significantly influences the survival time, against gravitational collapse, of the hyper-massive 
neutron star (HMNS) likely produced by the merger. It is therefore important that special 
attention is paid to use EOSs that are physically motivated, as done in 1431 within a 
conformally flat description of the fields and a simplified treatment of the hydrodynamics. 
Because we are here mostly concerned with drawing a first qualitative picture of the properties 
of the torus in a space of parameters that is as vast as computationally affordable, we have 
employed the commonly used "ideal-fluid" EOS, in which the pressure p is expressed as 
p = pe{T — 1), where p is the rest-mass density, e is the specific internal energy, and T is 
the adiabatic exponent. Such an EOS, while simple, provides a reasonable approximation 
and we expect that the use of realistic EOSs would not change the main results of this work. 
Furthermore, it was shown in l44l that the inspiral of equal-mass binaries of NSs described 
by realistic EOSs can be reproduced quite well by studying NSs with the same mass and radii 
but constructed as poly tropes with F = 2. 



Accurate evolutions of unequal-mass neutron-star binaries 



5 



As in Is), the gravitational-wave signal is extracted using two methods. The first method 
uses the Newman-Penrose formalism so that the gravitational-wave polarization amplitudes 
h+ and /ix are then related to ^'4 by simple time integrals 1451 

h+-ihx = 1'4 , (1) 
where the double overdot stands for the second-order time derivative and the curvature scalar 

*4 -Cap^5n"m^n''fh^ (2) 

is defined as a particular component of the Weyl curvature tensor, Cap-yS, projected onto a 
given null frame {l,n,m,Th} (see 1321 for details). The second and independent method 
is instead based on the measurements of the non-spherical gauge-invariant perturbations of a 
Schwarzschild BH (see refs. B6ll47ll48l for some applications of this method to Cartesian- 
coordinate grids). In particular, the gravitational-wave polarization amplitudes are in this case 
expressed in terms of gauge-invariant metric perturbations B9l 

- i/^x - ^ E - i W QL(^')rft') -2^^" , (3) 

where _.-2Y^'^ are the s = —2 spin-weighted spherical harmonics and Q^,„, the (gauge- 
invariant) odd-parity (or axial) current multipoles and even-parity (or polar) mass multipoles 
of the perturbed metric, respectively. In practice, these multipoles are computed on a set of 
2-spheres of fixed coordinate radius rjso = 200 {i.e. ^ 300 km). 

The two wave-extraction methods yield results whose differences are smaller than 1%, 
and hereafter we will concentrate only on the one using the gauge-invariant perturbations 
as it reduces the number of integration constants to be determined when computing the 
gravitational-wave strain. 



2.2. Adaptive Mesh Refinements and Grid setup 

The grid hierarchy is handled by the Carpet mesh refinement driver f50]. It implements 
vertex-centered mesh refinement, also known as the box-in-box method, and allows for 
regridding during the calculation as well as multiple grid centers. With mesh refinement, a 
small number of grids with different resolutions, called refinement levels, overlay each other, 
and are nested in a way that the coarsest grid has the largest extent and the finest grid the 
smallest extent. While the refined grids in the interior allow for an increased resolution where 
it is most desired, the outer boundary can at the same time be kept at a large distance. 

The timestep on each grid is set by the Courant condition (expressed in terms of the speed 
of hght) and so by the spatial grid resolution for that level; the typical Courant coefficient 
is set to be 0.35. The time evolution is carried out using 4th-order-accurate Runge-Kutta 
integration algorithm. Boundary data for finer grids are calculated with spatial prolongation 
operators employing 3rd-order polynomials and with prolongation in time employing 2nd- 
order polynomials. The latter allows a significant memory saving, requiring only three 
timelevels to be stored, with little loss of accuracy due to the long dynamical timescale relative 
to the typical grid timestep. 

For the inspiral phase of the system of binary NSs, two grid centers {rc.i : i = 1, 2} are 
defined, with one grid center located at the grid point where the rest-mass density reaches its 
maximum p,„ax = niax{p) and the other grid center located at the 7r-symmetric point (i.e. the 
grid point correspondent to pmax and rotated by 180 degrees around the z-axis). The grid 
hierarchy is composed of six refinement levels and a 2 : 1 refinement factor for successive 
levels. Once the condition pmax = niax(pi„ax,i) > 1-2 Pmax.initiai is satisfied, which is 
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known from experience to occur during the merger phase and well before collapse, the grid 
hierarchy is reduced to a single grid center fixed at the origin of the grid. At the initial time, 
the finest grids cover each star completely. Later, during the merger phase, matter outflows 
cross the boundary to the second finest grid and subsequently to the other coarser refinement 
levels. The grid resolution varies from Ag ~ 0.15 {i.e. ~ 221m) for the finest level to 
Ai = 4.8 {i.e. ^7.1 km) for the coarsest level, whose outer boundary is at 240 in our units 
{i.e. ~ 360 km). Initially, the number of grid points across the linear dimension of a star is of 
the order of 100. The torus surrounding the BH after collapse is usually not contained within 
the finest grid, but its high-density region is covered by the second finest grid with resolution 
As = 0.3. 

The whole grid is set up to be symmetric with respect to the {x, y) plane for both unequal- 
mass binaries and equal-mass binaries. The boundary conditions are chosen to be radiative 
for the metric, in order to prevent gravitational waves from scattering back into the grid, and 
static for the hydrodynamical variables. Note that the above setup is identical to that adopted 
inlH. 

Table 1. Properties of the binary NS initial data. From left to right the columns show: the 
name of the model (assembled from its rounded total baiyonic mass preceded by the letter 
M and its mass ratio preceded by the letter q), the total baryonic mass A/tot of the system, 
the total ADM mass A/^j^^, of the system, the ratio of the baryonic masses of the two stars 
q = M2/A/1, the ratio of the ADM masses of the two stars, the total angular momentum 
J, the initial orbital frequency Vorb^ the initial maximum rest-mass density pmax, the mean 
radius f,; of each star, the axis ratio Ai of each star, the individual ADM mass M°° of each 
star as considered in isolation at infinity, and the compactness C°° of each star as considered 
in isolation at infinity. The mean radius is defined as fi = (r^ + r-{ + rj_ + ''pol)/4, where 
and are the radii of the star parallel to the line connecting the stars, rj_ is the radius 
in the equatorial plane perpendicular to that line, and rpoi is the radius perpendicular to the 
equatorial plane. The axis ratio is defined as the ratio between the mean radius parallel to the 
line connecting the stars, and the mean radius in the plane perpendicular to that line, namely 
Ai = {r± + rpai)/{r\~ + r-\). All values except pmax are provided by the output of the 
LORENE code, and the accuracy of Mtot and J is the one at which the Whisky code is able 
to reproduce them for the present setup. 
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13.1, 12.1 
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. 4q0 . 
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8.36 
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2.3. Initial data 

We use quasi-equihbrium initial data generated with the multi-domain spectral-method code 
LORENE developed at the Observatoire de Paris-Meudon f5T\. For more information on the 
code and its methods, the reader is referred to the LORENE web pages [52]. In particular, 
because the binaries are not expected to be corotating |f53l , we use irrotational configurations, 
defined as having vanishing vorticity, and obtained under the additional assumption of a 
conformally flat spacetime metric |51il . 

Some of the models investigated in this paper are publicly available on servers of 
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the Meudon group ll52ll . Others have been created by us specifically for the unequal- 
mass simulations presented here. The models of the lowest mass ratios have been kindly 
provided by Dorota Gondek-Rosiiiska). The EOS assumed for the initial data is in all cases 
the polytropic EOS p = K with an adiabatic index F = 2 and a polytropic coefficient 
K = 0.0332 pnuc c^/j^nuc — 123.6 (in units where c = G = Mq = 1), where /9„uc and n^^c 
refer to the nuclear rest-mass and number densities, respectively. For this particular EOS, 
the allowed maximum baryonic mass for an individual stable NS is ^ 2.00 M0. The initial 
coordinate separation of the stellar centers in all cases is c? = 45 km. 

The models used as initial data include both equal-mass models and most importantly 
unequal-mass models. As mentioned in the Introduction, we here concentrate on the dynamics 
of the massive tori resulting from the merger. The formation of the torus is enhanced for 
smaller mass ratios. Although the observational evidence is not very firm ifTSl 1541 there 
is also no theoretical reason ruling out their possible existence ifTSl [TtII . A full list of all 
considered models together with a selection of physical quantities defining them, e.g. baryon 
and ADM mass, orbital frequency and initial angular momentum, etc., is given in table [T] 
To distinguish simply the different binaries we adopt the following naming convention: any 
initial data binary is indicated as M%q#, with % being replaced by the rounded total baryonic 
mass Mtot of the binary neutron-star system and # by the mass ratio q. As an example, 
M3 . 4 qO . 8 is the binary with total baryonic mass Mtot — 3.4 Mq and mass ratio q = 0.80. 

3. Dynamics of the coalescence and merger 

3.1. General dynamics 

In a previous work 13], we have investigated the dynamics of the coalescence and merger 
of equal-mass binary NSs for models with total baryonic mass A/tot — 2.912 Mq and 
-Mtot = 3.250 Mq. It was found that for any of the two EOSs considered, binaries with 
(initial) total baryonic mass below a certain limit do not collapse promptly to a BH but rather 
yield an oscillating HMNS, which undergoes delayed collapse to a BH. Independently of 
the mass ratio, all of the binaries under consideration here have masses higher than those 
considered in fSj and all collapse promptly never leading to a HMNS even if the EOS used 
here is a non-isentropic one (see discussion in |fj| on the different qualitative behaviour 
between an isentropic and a non-isentropic EOS). This absence of a HMNS, however, is the 
direct consequence of the chosen initial data rather than a feature of unequal-mass mergers 
and it has been here exploited simply to reduce the computational costs by shortening the time 
of the collapse to a BH. 

Figure [T] shows a selection of representative isodensity contours on the equatorial plane 
for the equal-mass binary M3 . 6ql .0 0. At the initial time, the stars are in their quasi- 
equilibrium configuration at a coordinate separation of 45 km. The binary progressively 
speeds up while inspiralling. After slightly more than two orbits have been completed (namely 
after about 5-6 ms), the stars merge and, about 2-3 ms later, an apparent horizon (which we 
search with the code of ll55l ) is found. The ideal-fluid EOS employed in the simulations allows 
for shock-heating and an increase of the specific internal energy e, as shown in Q; this, in 
turn, causes some matter to be ejected from the rotating central object and to propagate into 
the surrounding artificial atmosphere. The evolution of model M3 . 6ql . shows that matter 
is ejected in small amounts during the inspiral phase and in larger amounts during the merger 
phase, when the shocks are much stronger Therefore, while small spiral arms can certainly be 
observed in the outer regions during the merger phase (see the last two snapshots of figure[T]i, 
they do not have sufficient angular momentum to reach distances as large as in the unequal- 
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Figure 1. Isodensity contours for the M3 . 6ql . model on the {x, y) plane. The times when 
the frames have been taken are shown on top of the plots while the color-code for the rest-mass 
density is indicated to the right of each plot. Additionally, isodensity contours are shown for 
the values of p = 10^", lO^S lO^^ ;^oi2.5 ;^gi3 j^qI-'-^, 10^"*, lO"'^, lO^''^ g/cm'^ The 
third frame (at time t = 5.760 ms) shows the onset of the merger, the last two frames (at times 
t = 8.210 ms, t = 8.276 ms) show the behaviour of the system during the collapse to a BH. 



mass models (see discussion below). Instead, the spiral arms wind around the rapidly rotating 
central object formed by the two NS cores. Quantitative results regarding the BH spin and the 
mass and angular momentum of the remaining disk will be discussed in subsequent sections. 
To contrast the evolution of an unequal-mass binary, figure |2] shows the same selection 
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Figure 2. Isodensity contours for the M3 . 4q0 . 7 model on the {x, y) plane. The times when 
the frames have been taken are shown on top of the plots while the color-code for the rest-mass 
density is indicated to the right of each plot. Additionally, isodensity contours are shown for 
the values of p = 10^", lO^S lO^^ ;^oi2.5 ;^gi3 ^qI^-S, IQi"*, lO"'^, lO^''^ g/cm'^ The 
third frame (at time t = 4.104 ms) shows the onset of the merger, the last two frames (at times 
t = 6.620 ms, t = 7.414 ms) show the behaviour of the system during the collapse to a 
BH. Note that the computational domain is much larger than what is shown and extends to 
~ 360 km 



of isodensity contours on the equatorial plane as represented in figure [T] only now for 
the M3 . 4qO . 70 model, which has the smallest mass ratio considered in this work. The 
asymmetry of the binary system is already apparent at the initial time. The heavier star is 
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much more compact than its extended less massive companion, which is deformed already 
at the initial distance by tidal forces. In addition, the center of mass does not coincide with 
the point halfway between the centers of the stars, but it is shifted toward the more massive 
star. During the inspiral phase, the heavier and more compact star is only slightly affected 
by its companion, whereas the latter is decompressed rapidly while being accreted onto the 
heavier star This is visible in the three intermediate panels of figure |2] The tidal disruption 
of the lower-mass NS when it still retains a large fraction of its angular momentum results in 
an extended tidal tail, which, unlike what happens in the equal-mass case, transfers angular 
momentum outwards in a much more efficient way. This leads to the formation of large spiral 
arms extending well beyond the domain shown in figure |2] and ultimately to a more rapid 
ejection of matter. Gravitationally bound matter travelling along the spiral arms away from the 
central object will form a more massive accretion torus around the central BH than that formed 
in the case of an equal-mass, symmetric binary system. It should be noted that, although the 
rest-mass density of the matter in these spiral arms is much smaller than the central one, it has 
nevertheless densities p > 10^" g/cm and thus well in a general-relativistic regime. 

3.2. Properties of the black hole 

As mentioned above, because of the large initial mass of the system and irrespective of the 
mass ratio, the merged object rapidly collapses to a BH. Its mass and angular momentum have 
been computed making use of the dynamical-horizon formalism 15611571 . which provides a 
simple and accurate measure of the BH properties also when this is subject to the inflow 
of mass and angular momentum fJSl. In the case of the equal-mass binary, because the 
disk resulting from the merger has comparatively small mass, the BH settles rapidly to 
an approximately stationary configuration, and the mass and spin of the BH measured at 
formation, i.e. M = 2.56 Mq and a = J /Al"^ = 0.745, respectively, do not vary significantly 
throughout the subsequent evolution of the system. On the other hand, when considering 
unequal-mass binaries, the mass and the spin of the BH show, on the timescale of the 
simulations, a variation in time of ~ 5% and ~ 2%, respectively, because of the continued 
and intense accretion of both mass and angular momentum. Table |2] shows the corresponding 
parameters for all models at the final time of the evolution, which is not the same for the 
different binaries considered. 

We note that finding and tracking the apparent horizon in the case of binaries with small 
mass ratio is far from being simple since the asymmetry in the merger dynamics leads to 
a noticeable motion of the "center-of-mass" of the system. Hence, the location of the trial 
surface for the apparent horizon cannot be simply associated to a pre-existing black hole (as 
in the case of BH binaries ll32l ) or to a pre-determined coordinate location (as in the case of 
the collapse of a rotating star (351). The end-result of this complication is that the apparent 
horizon could not be tracked successfully in all the models under consideration. This was the 
case of models M3 . 4q0 . 80 and M3 . 4q0 . 7 0, for which it was not possible to measure the 
mass and spin of the corresponding BH. Furthermore, in the case of the binary M3 . 5q0 . 7 5 
the measurements were not made with the dynamical-horizon formalism but rather by using 
the ratio of the polar to equatorial circumference of the apparent horizon as discussed in detail 
in 1351 . Cross-checking the two measures (i.e. apparent-horizon distortion and dynamical- 
horizon formalism) in the cases where both are possible shows that they are equally reliable 
(see also the extended discussion in l35l ). Overall, the data available suggests the existence 
of a local maximum of a for q ^ 0.9, but more data is clearly necessary to confirm this. 

Interestingly, when inspecting carefully the apparent horizon in the low-q model 
M3 . 5qO . 75 it is possible to appreciate that its appearance precedes the time when the 
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two stellar cores merge and is in contrast with what happens with models with high q. By 
comparison, we believe the same happens also for the binary M3 . 4q0 . 7 0, although in this 
case we were not able to detect an apparent horizon. Of course these considerations have little 
physical importance as the interior of the apparent horizon is causally disconnected with what 
is astrophysically observable; nevertheless this result provides another interesting example of 
the rich phenomenology relative to the appearance and dynamics of trapped surfaces (see, for 
instance, the discussion in section 4 of Ii58j| ). 

Table 2. Columns 2 — 5 report the properties of the final BH, i.e. mass, angular momentum, 
spin parameter, and kick velocity, while columns 6 — 7 report the measured torus masses Aftor 
and those inferred from relation j6j, Mtor ■ Also shown in columns 8 and 9, respectively, are 
the numerical error etor = |1 — (MtoiOfjR / (^^tor)^!? I computed by comparing different 
resolutions (medium, i.e. Ag = 0.19, and high, i.e. Ag = 0.15) for each model and the 
relative error egt = |Aftor — Aftor|/A/tor of the phenomenological expression for the mass 
of the torus with respect to the numerical data. Clearly, the binaries with high mass ratio are 
not well described by relation (6) even though their numerical eiTor is not very lai'ge. 



Model M J a = J/KP iJkick Mtor |A'/toi-| etor tflt 

(km/s) (Mo) (Mo) 



M3 . 


. 6ql . 


.00 


2.56 


4.90 


0.745 


0.28 


0.0010 


0.021 


28% 


> 100% 


M3 . 


. 7q0 . 


. 94 


2.64 


5.18 


0.743 


121.95 


0.0100 


0.048 


12% 


> 100% 


M3 . 


. 4q0 . 


. 91 


2.99 


7.29 


0.815 


59.33 


0.0994 


0.103 


0.8% 


89.6% 


MS . 


, 4q0 . 


.80 








56.22 


0.2088 


0.193 


1.5% 


7.4% 


M3 . 


. 5q0 . 


. 75 


3.00^ 


7.13''' 


0.792t 


18.05 


0.0802 


0.173 


2.5% 


8.1% 


M3 . 


. 4q0 . 


.70 








15.82 


0.2116 


0.202 


2.4% 


4.6% 



' We could not compute the dynamical horizon for this model, so the reported values are calculated 
from the apparent horizon, with the method employed in Sections VA and VBl of ref. 1351 . 



Finally, also reported in table |2] is the recoil velocity imparted to the BH at the end of 
the inspiral and computed using the gravitational-wave emission as discussed in 1591 [32l for 
binary BHs. We recall, in fact, that together with energy and angular momentum, gravitational 
radiation also carries away linear momentum. If the binary system has a degree of asymmetry 
(either in the mass or in the spin) then the trajectories of the two bodies will be (slightly) 
different {e.g. with the smaller body moving more rapidly and, hence, being more efficient 
in beaming its emission) and the momentum loss in one direction will not be balanced by an 
equal loss in the diametrically opposite direction. This effect is well-known in binary BHs, 
where the recoils from quasi-circular inspirals can be as larger as ~ 4000 km (see 1601 for 
a recent review), but it has never been reported before for binary NSs. The recoil velocities 
reported in table |2] are clearly much smaller than those measured for binary BHs. However, 
they could still yield astrophysically interesting results being comparable or larger than the 
escape velocity from the core of a globular cluster that is Vcsc ^ 50 km ||6T1. Furthermore, 
and possibly surprisingly, the values reported here for irrotational binaries which have very 
little initial spin, are not much smaller than those computed for non-spinning binary BHs (see, 
e-g- , l62ll for a recent update) and have a local maximum for q ^ 0.9. 
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Figure 3. Isodensity contours for the binaries M3 . 6ql . 00 (left panels) and M3 . 4q0 . 70 
(light panels) showing the morphology of the toii at the onset of the QSA on the {x, y) 
plane (upper rows) and on the {x, z) plane (lower rows). Note that the disks in the 
two panels have very different lengthscales, with the one for M3 . 4q0 . 70 being about 3 
times larger than that for M3 . 6ql . 00. The colormap used here is different from the one 
in figures [T] and |2] Additionally, isodensity contours are shown for the values of p = 
10^", 10", 10^2, 10^=* g/cm^*. 



4. Torus Formation and Properties 

4.1. General Dynamics 

In figures [3] and |4] we show color-coded contours of the rest-mass density for models 
M3 . 6ql . 00 (left panels) and M3 . 4q0 . 70 (right panels), either in the (a;,y) plane (upper 
rows) and in the [x, z) plane (lower rows). The snapshots in figure[3] in particular, correspond 
to the time t ^ \Q ms when the systems enter the regime of quasi-stationary accretion (QSA, 
see below for definition), shortly after the formation of the BH, while those in figure |4] refer 
to the final time of the evolution, t ^ 2 \ ms. These figures allow for a closer view of the 
morphological features of the disks, in particular, their spatial dimensions and thickness, and 
are a natural continuation of the dynamics already shown in figures[T]and|2] although they use 
a different colormap that has been tuned to yield a better contrast in the density profiles. 

The large morphological differences between these two extreme models are clearly 
visible in these figures. The equal-mass model produces a highly symmetric, geometrically 
thin disk, similar to the ones already observed for other equal-mass initial data in |f3!|. The 
unequal-mass model, on the other hand, at the time of the onset of the regime of QSA is 
characterized by the presence of a large spiral arm, which has not yet been accumulated onto 
the central disk surrounding the formed BH. The asymmetry in the distribution of matter 
at this stage is also apparent from the color map of the rest-mass density. Only at the end 
of the simulation the disk of the unequal-mass binary acquires a more axisymmetric shape. 
The diameters of the disks and their heights perpendicular to the horizontal plane differ in a 
significant way between the two models. More specifically, at the end of the evolution and 
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Figure 4. The same as figure[3]but showing the tori at the end of the simulation. 



using the p = 10^° g/cm^ isodensity contour as the reference value below which material is 
not considered part of the disk, our simulations yield disk diameters of ^ 50 km for model 
M3 . 6ql . and ^150 km for model 1X13 . 4 qO . 7 . The corresponding vertical scale is ~ 5 
km and ~ 35 km, respectively^ Taking into account all the models of our sample, we find that 
both scales increase as the mass ratio decreases. Even more worth noticing is the fact that, in 
the cases considered, while the tori differ in size by about a factor ~ 3, they differ by a factor 
200 in mass, while having comparable mean rest-mass densities (see further discussion in 
sections |431 and I4l8b. 



4.2. Rest-mass Evolution 

In order to establish how the asymmetry in the mass of the two NSs in the binary leads to tori 
with different masses we show in figure|5]the evolution of the total rest mass, defined as 

Aftot= j pW^£x = j D^d'x, (4) 

V V 

normaUzed to its initial value and for the different models. In this equation W = au* is 
the Lorentz factor, a being the lapse function, and 7 is the determinant of the spatial metric. 
All the curves in figure [5] have been shifted in time to coincide at tcow, which represents 
the (collapse) time at which a rapid decrease of the total rest-mass takes place following the 
formation of a BH. Note that, in practice, the collapse time is different for all models, ranging 
from around 6 ms for model M3.4q0.70to around 11 ms for model M3 . 4 qO . 8 . 

X Of course it should be noted that the spatial dimensions reported here depend on the cut-off chosen for the rest- 
mass density. Using cut-offs smaller than p = 10^" g/cm'' would lead to considerably larger estimates for the sizes 
of the tori. 
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Figure 5. Evolution of the total rest masses A/tot normalized to their initial values for all the 
models considered. The order of magnitude of the mass fraction in the accretion torus can 
be read off the logarithmic mass scale on the vertical axis. The curves refeiTing to different 
models have been shifted in time to coincide at tcoih which represents the time when the very 
rapid decrease of the total rest mass takes place. Note that this time is not physically relevant 
(the apparent horizon is usually found earlier than tcoll) and simply coiTesponds to when the 
very large amount of rest-mass accumulated in a very few cells is numerically dissipated. 



Figure ID shows that all models conserve the baryonic mass almost perfectly (i.e. with 
losses of < 10~^) up until the formation of the apparent horizon, after which most of the rest 
mass disappears in the singularity. One obvious result which can be deduced from figure|5]is 
that the mass of the resulting accretion disk is larger for smaller values of q. However, this 
trend is not entirely monotone in the figure as it is also influenced by the initial total baryonic 
mass of the binary. The particular values of the tori masses computed for all models are 
reported in table|2] While the equal-mass model produces a disk of barely 10^"^ A/0, models 
M3 . 4q0 . 8 and M3 . 4q0 . 7 produce significantly more massive tori with masses of about 
0.2 Mq. a more detailed discussion of the mass in the tori will be made in section |4~8] 

As the apparent horizon is formed, a substantial part of the rest-mass is still outside it, 
although it will accrete rapidly onto the BH. This makes the mere definition of what is the 
torus and its mass rather arbitrary and we decide therefore to define the torus mass Aftoi as 
the total rest mass outside the apparent horizon when the disk enters a regime of quasi-steady 
accretion (QSA), a regime which is found in all models investigated. More specifically, we 
compute the accretion rate as 

7 p 

Mtot^j^ pW^d'x (5) 

V 

and define the onset of the QSA as the point in time when the condition A/tot /A/tot < 
lQ~^{Gc~^ Mq)~'^ is satisfied for the first time. In other words, we define the onset of 
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t (ms) I (ins) 

Figure 6. Evolution of the total baryonic mass J\/tot (upper rows) and of the accretion rate 
A/tot (lower rows) in the regiine of QSA for the representative models M3 . 6ql . 00 (left 
panel) and M3 . 4q0 . 70 (right panel). Indicated with a vertical dashed line is the onset of 
the QSA. Note that both A/tot and Aftot differ by almost two orders of magnitude when 
coinparing equal and unequal-mass binaries. 

the QSA as the time when the accretion has stabiHzed and the matter is moving on essentially 
circular orbits. This definition is again somewhat arbitrary, but has the advantage of allowing 
for a systematic comparison of the differences in the properties of the accretion tori produced 
by the several models considered in this work. 

Figure |6] shows the evolution of the total rest mass A/tot (upper rows) and the mass 
accretion rate A/tot (lower rows) in the regime of QSA for the two extreme models of our 
sample, M3 . 6ql . (left panel) and M3 . 4q0 . 70 (right panel). Also indicated with a vertical 
dashed line is the onset of the QSA and it should be noted that both A/tot and A/tot differ by 
almost two orders of magnitude when comparing equal and unequal-mass binaries. An aspect 
of the evolution of the accretion rate which is quite evident in figure |6] is the sharp difference 
between equal- and unequal-mass binaries. The equal-mass case, in fact, shows an accretion 
rate (and indeed the whole evolution of the torus) that is subject to quasi-periodic oscillations 
as the torus moves in and out at about the radial epicyclic frequency. The mass flux of the 
unequal-mass model, on the other hand, is rather constant in time and this reflects a very 
different distribution of angular momentum in the tori. Both of these aspects will be further 
discussed in the following sections. 

4.3. Density evolution 

Once the BH is formed, an effective gravitational potential well builds up in which the 
torus undergoes radial oscillations. In the case of an equal-mass binary, the well is 
essentially axisymmetric and the dynamics of oscillating relativistic tori in equilibrium and 
in axisymmetry have been analyzed extensively in a series of papers 1631 l28l [30l l26l in the 
test-fluid approximation (where the self-gravity of the disk is neglected), with and without 
magnetic fields, and for the cases of Schwarzschild and Kerr BHs. These papers have shown 
that, upon the introduction of perturbations in the tori, a long-term oscillatory behavior is 
found, lasting for tens of orbital periods. These oscillations correspond to axisymmetric 
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Figure 7. Evolution of the maximum of tlie rest-mass density pmax, normalized to its initial 
value for the representative models M3 . 6ql . 00 (left panel) and M3 . 4q0 . 70 (right panel). 
The rapid drops take place well after an apparent horizon has been formed and are caused by 
the numerical methods which are no longer able to resolve the very large gradients in the very 
central grid cells. The two insets provide a magnified view of the evolution of the density in the 
torus and help to contrast the periodic accretion produced in the case of equal-mass binaries 
and the QSA for the unequal-mass binaries. 



p-mode oscillations whose lowest-order eigenfrequencies appear in the harmonic sequence 
2:3. This harmonic sequence is present with a variance of ^ 10% for tori with a 
constant distribution of specific angular momentum and of ^ 20% for tori with a power- 
law distribution of specific angular momentum. More recently, those studies have been 
extended by 11641 . where systems formed by a BH (in the puncture framework) surrounded by 
(marginally stable) self-gravitating disks have been evolved in axisymmetry. Even in this case, 
the ratio of the fundamental oscillatory mode and the first overtone also shows approximately 
the 2:3 harmonic relation found in earlier works ll28l [30l l26l . 

The dynamics of the BH-torus system produced by the merger of model M3 . 6ql . is 
considerably more complicated than that considered in the test-fluid studies, for which initial 
configurations in stable equilibrium could be found. However, despite the fact that these 
systems are studied ab-initio as the end-products of highly dynamical events, it is remarkable 
that so much of the phenomenology studied and reported in [|63..28..30. ,261 continues to apply 
also here. Unfortunately, although the simulation extends to ~ 26 ms, the timeseries is much 
too short to provide a firm evidence of the presence of the 2:3 harmonic relation, although the 
spectral analysis of the data indicates that excess power is present at such frequencies. 

To provide additional evidence that the harmonic behaviour is not just in the accretion 
rate, figure Q shows the evolution of the maximum of the rest-mass density p,„ax, normalized 
to the corresponding initial value, for the two extreme models of our sample, M3 . 6ql . 00 
andM3 . 4q0 . 70. The equal-mass model M3 . 6ql . 00 (represented in the left panel of figure 
shows the most regular and pronounced oscillatory behavior, as was already evident in the 
time evolution of the total rest mass and accretion rate in figure|6] The two insets in this figure 
magnify these features in the QSA regime and, in the case of M3 . 6ql . 0, they highlight the 
presence of both maxima and minima corresponding to configurations when the torus reaches 
the point of closest approach to the BH (periapsis or pericenter) and of farthest excursion 
(apoapsis or apocenter), respectively. A similar trend can be hinted also for the M3 . 4q0 . 7 
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model on the right panel of Fig [T] although the quality of the oscillations is smaller in this 
case, most likely because in this case the enhanced tidal disruption during the merger phase 
leads to a more complex dynamics. Interestingly, such oscillations seem to become more 
regular during the final stages of the evolution, i.e. for t > 17 ms, as the torus reaches a more 
axisymmetric configuration. 

A novel technique to analyze the evolution of the tori and to gain some insight on their 
dynamics is that offered by spacetime diagrams for observers comoving with the black hole. 
This is shown in figure [8] which reports the evolution of the color-coded rest-mass density 
embedded in a spacetime diagram with the (x — x^^^ ) coordinate on the horizontal axis, where 
x^^ is the position of the apparent horizon, and the coordinate time t on the vertical axis. The 
color code is indicated to the right of each plot and isodensity contours are shown for the 
values of p = 10^°, 10^\ 10^^, 10^^ g/cm^. Note that while these values are the same for all 
the panels, the spatial dimensions vary considerably. For each model, the dotted horizontal 
line marks the onset of the regime of QSA. 

By comparing the spacetime diagrams for all models it is evident that only the equal-mass 
model M3 . 6ql . 00 shows a global oscillatory movement with respect to the location of the 
BH horizon. The movement is indeed global as all the isodensity contours plotted oscillate 
simultaneously and the maximum and minimum radial extensions reached by the disk (as 
signalled by the location of the 10^° g/cm'' contour) are ~' 25 km and ^ 15 km, respectively. 
It is these oscillations that produce the periodic increase in the maximum rest-mass density 
reported in figureQand it is easy to appreciate that in this case the average density in the disk 
is less than about 10^^ g/cm"^. 

Scrolling through the different panels in figure |8] it is possible to appreciate that the 
dynamics of the torus is strongly influenced by the mass ratio. More specifically, models 
M3.4q0.80,M3.5q0.75,M3.4q0.70, show very rapid expansions corresponding to the 
ejection of the large spiral waves discussed in the previous sections. As we will comment 
later on, most of this matter is still bound but it nevertheless reaches distances which are 
several hundreds of km away from the BH, leading to tori that have spatial dimensions as 
large as ~ 80 km. Furthermore, noticeably higher average rest-mass densities are reached in 
the three low-q models. As a result, while tidal disruption sweeps away a large fraction of the 
external layers of the less massive star in the binary, the corresponding tori are still able to 
retain the inner and denser regions; this is particularly the case for models M3 . 4q0 . 80 and 
M3 . 4q0 . 7 0, where the tori reach maximum densities as high as ^ 10^"* g/cm^. 

We finally note that, because the matter in the spiral arms is bound, it will eventually fall 
back onto the tori, where it may lead to enhanced accretion and consequently to a new and 
delayed gamma-ray emission as the one recently observed in Il65ll . Determining the dynamics 
of the fall-back material is therefore of great importance astrophysically and the focus of our 
attention in future work. 

4.4. Dynamical Instabilities 

As mentioned in the Introduction, current models of GRBs assume that the central engine is a 
system consisting of a BH and a thick disk, either formed at the late stages of the coalescence 
of two NSs or after the gravitational collapse of a massive star. The energy supply comes from 
the energy released by the accretion of disk material on to the BH and from the rotational 
energy of the BH itself, which can be extracted, for instance, via the Blandford-Znajek 
mechanism f66l. This vast amount of energy (of the order of lO^'^-lO^* erg, depending 
on the mass of the disk and on the BH rotation and mass) is sufficient to power a GRB if 
the energy released can be converted into 7-rays with an efficiency of about a few percent. 




Figure 8. Evolution of the rest-mass density p along the positive x axis in a frame comoving 
with the BH. The panels show the color-coded rest-mass density embedded in a spacetime 
diagram with the [x — Xj^^) coordinate on the horizontal axis, being x^jj the position 
of the apparent horizon, and the coordinate time t on the vertical axis. The color code is 
indicated to the right of each plot. Additionally, isodensity contours are shown for the values 
of p = 10^°, 10^\ 10^^, 10^3 g/cm^. For models M3 . 4q0 . 80 and M3 . 4q0 .70, where 
the horizon could not be tracked, Xj^-^ represents a guess for the border of the horizon. For 
each model, the dotted horizontal line marks the onset of the regime of QSA. 
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This scenario requires a stable enough system to survive for a few seconds. In particular, the 
internal-shock model 1671 implies that the duration of the energy release by the source has 
a duration comparable with the observed duration of the GRB. Any instability which might 
disrupt the system on shorter timescales, such as the so-called runaway instability |68 1, could 
pose a severe problem for the accepted GRB models. The runaway instability was first pointed 
out in 1201 and operates as follows: If the torus is initially filling its Roche lobe, transfer of 
mass onto the BH is possible through the cusp located at the Li Lagrange point. As a result 
of accretion, the mass of the BH increases, thus leading to a change in the gravitational field 
of the system and ultimately to a change in the position of the cusp. This can move either 
inwards (towards the BH) or outwards (away from the BH), and when the latter happens it 
leads to an increase in the mass transfer and hence to the runaway accretion of the torus on a 
timescale of a few milliseconds. 

The runaway instability has been investigated under different assumptions and 
approximations (see l22l l24l and references therein). Early simplified studies based on 
stationary models showed that, on one hand, the self-gravity of the disk favours the instability, 
and, on the other hand, there are also parameters which may help to stabilize the disk, such 
as the rotation of the BH and the radial distribution of specific angular momentum. The 
first time-dependent, general relativistic hydrodynamical axisymmetric simulations of the 
runaway instability of tori around BHs were performed by 122 , 69l l25ll23l . who treated the 
dynamics of the gravitational field in an approximate way and neglected the self-gravity of 
the torus. Overall, ll22l l25l |23]| found that tori with constant distribution of specific angular 
momentum were unstable while non-constant (power-law) angular momentum disks were 
stable. More recently, in ll24ll the first simulations in full general relativity of marginally- 
stable self-gravitating tori in axisymmetry were performed with the purpose of evaluating the 
influence of the torus self-gravity on the runaway instability. The results of l24l indicate that 
the tori are indeed stable irrespective of the angular momentum distribution. It is therefore 
interesting that the results presented in figure |8] which are not restricted to axisymmetry 
but are however constrained to much shorter timescales, reach the same conclusion: Self- 
gravitating tori around BHs, as those produced by the merger of binary NSs, are stable at least 
on the dynamical timescales investigated here. Additional considerations on the stability of 
the tori are presented in the following section. 

4.5. Specific Angular-Momentum Evolution 

Besides the rest-mass density, another quantity whose evolution is useful to understand the 
dynamics of the tori is the the specific angular momentum. This quantity plays an important 
role in defining the dynamics of point particles around black holes and in defining the 
equilibrium of non-self gravitating tori around black holes fTOl. As mentioned above, we 
define the specific angular momentum as ^ = — U0 /uf. We also note that a similar but distinct 
definition of the specific angular momentum was used in |1|, namely j ~ hu^. The two 
definitions have the same Newtonian limit of jNcwt = ^Nowt = ilr^, being the angular 
velocity. However, it is important to stress that only the definition used here yields the correct 
zero radial epicychc frequency for tori with constant specific angular momentum [see eqs. (43) 
and (45) of l28ll1. 

Figure |9] shows the evolution of the specific angular momentum, where the different 
panels show the color-coded specific angular momentum for an observer comoving with 
the BH. The color code is indicated to the right of each plot and in addition the same 
isodensity contours reported in figure |8] are shown here to aid to follow the dynamics of 
the matter The most striking feature to note when scrolling through the different panels in 
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Figure 9. The same spacetime diagrams as in figure [8] but for the evolution of the specific 
angular momentum £ = —u^/ut- Note that the isocontours in this case refer to the rest-mass 
density and are the same as in figure[8] 



figure |9] is that the radial distribution changes radically but systematically when going from 
the equal-mass binary M3 . 6ql . over to the most extreme unequal-mass binary considered 
M3.4q0.70. In particular, while the specific angular momentum is decreasing outwards 
in models M3 . 6ql .00, M3 . 7qO . 94, and M3 . 4qO .91, it is Keplerian and increasing 
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Figure 10. Profiles along the x-axis of the specific angular momentum of the tori produced 
by the binaries M3 . 6ql . (blue lines extending to < 20 km) and M3 . 4q0 . 7 (red lines 
extending up to > 70 km). The profiles are computed in a frame comoving with the BH 
and for densities p > 10^'' g/cm'^. Different line types refer either to the onset of the QSA 
(i.e. t ~ 10 ms, thin solid lines) or to the end of the simulation (i.e. t ~ 22 ms, thick dashed 
lines). Note the markedly different behaviour and that the specific angular momentum for the 
unequal-mass case increases outwards. 



outwards as ^ x^/^ for the remaining models (see also the discussion in the following 
section). Furthermore, the spacetime plots show that the matter located in the outer regions 
of the disks acquires the largest values of the specific angular momentum. This is particularly 
visible in the early evolution of model M3 . 4q0 . 80, in which a large spiral arm develops 
extending beyond the computational boundary, and also in the late evolution of model 
M3 . 4q0 . 7 0, when the corresponding disk reaches the largest radial extension {cf. right 
panel of figures [3] and |4]i. Broadly speaking, our simulations show that, in agreement with 
the results of [Tl, the smaller the value of q, the more the angular momentum is transported 
outwards by a torque from the non-axisymmetric object that forms after the merger 

To better highlight the different behaviour of £ for different mass ratios we show in 
figure[TO]the profiles along the x-axis for the tori produced by the binaries M3 . 6ql . (blue 
lines extending to < 20km) and M3 . 4q0 . 70 (red lines extending up to > 70km). The 
profiles are computed in a frame comoving with the BH and for densities p > 10^" g/cm^. 
Different line types refer either to the onset of the QSA (i.e. t ~ 10 ms, thin solid lines) or 
to the end of the simulation {i.e. t ~ 22 ms, thick dashed lines). Quite clearly, the specific 
angular momentum decreases outward at all times for the equal-mass binary, while it increases 
outward for the unequal-mass one (although it was initially decreasing at the innermost parts). 
At this point it is worth remarking that Rayleigh's criterion against axisymmetric perturbations 
of rotating fluids requires that d£/dx > for a dynamical stability |71|. While this criterion 
is clearly satisfied by model M3 . 4q0 .70, it is equally-clearly violated by M3 . 6ql . 00, 
which is nevertheless stable. We believe this difference is due to the fact that Rayleigh's 
criterion assumes that the motion is stationary and purely azimuthal. While this is essentially 
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the case for the unequal-mass binary, which does not show clear evidence of epicyclic 
oscillations, it does not hold true for the unequal-mass binary, which shows instead large radial 
epicyclic oscillations. The nonlinear stability of M3 . 6ql . 0, but also of M3 . 7qO . 94 and 
M3 . 4q0 . 91, seems therefore to indicate that Rayleigh's criterion can and should be extended 
to account for fluids which are subject to large radial excursions. 



M3.6q1.00 



M3.7q0.94 




Figure 11. The same spacetime diagrams as in figure|9]but for the evolution of the angular 
velocity Q. Note that the isocontours in this case refer to the rest-mass density and are the 
same as in figure[8] 
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Figure 12. The same as in figure \W\ but for the angular velocity. Shown as reference with 
a dotted line is the Keplerian angular velocity Qj^^p, which matches very well the outer 
parts of the torus from the unequal-mass binary. Shown instead with a long-dashed line is an 
exponentially decaying profile, which instead reproduces well the profile for the equal-mass 
binary. 



4.6. Angular-velocity Evolution 

In analogy with figures [8] and |9l figure [TT| shows the spacetime diagram for the evolution of 
the angular velocity ^l = u'^/u* for all models of our sample. It is straightforward to notice 
that for all models the angular velocity decreases with the radial distance from the apparent 
horizon. While this is qualitatively in agreement with the results of H], it is worth noting that 
the radial fall-off is very different as the mass ratio is varied among the different binaries. This 
is shown in figure [12] which reports the profiles of fl along the x-axis for the tori produced 
by the binaries M3 . 6ql . 00 (blue lines extending to < 20km) and M3 . 4q0 . 70 (red lines 
extending up to > 70 km). As before, the profiles are computed in a frame comoving with 
the BH and for densities p > 10^°g/cm^ and different line types refer either to the onset 
of the QSA [i.e. t ^ 10 ms, thin solid lines) or to the end of the simulation [i.e. t ^ 22 ms, 
thick dashed lines). It is then clear that while the equal-mass binary has an exponentially 
decaying profile (cf. long-dashed line), i.e. ft cx cxp[— fc(a; — .t^^j )] ~ cxp[— 0.07(a; — x^^^)], 
which does not change significantly with time, the unequal-mass binary reaches at the end of 
the simulation a profile which is, especially in the outer parts, essentially Keplerian, i.e. with 
^Kcp ^ x""^^^ (cf. dotted line). This feature, which is also shared by the other low-q binaries, 
explains the scaling of the specific angular momentum as ^ ~ x^^^ and provides firm evidence 
that the tori produced in this case will be dynamically stable. 

4.7. Matter Ejection 

As a final but nevertheless important aspect of the formation and evolution of the tori, we 
consider whether or not a part of the rest-mass of the system is ejected during the merger 
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and the subsequent evolution. To determine whether a fluid particle is bound or unbound we 
use the covariant time component of the 4-velocity Ut and recall that, in an axisymmetric and 
stationary spacetime, the value of ut for a particle moving along a geodesic is conserved. If the 
particle is unbound, it moves outwards and —ut ~ > 1 at infinity, where a = 1, = 0. 
The local condition > — 1 thus provides a necessary although not sufficient, condition for 
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a fluid element to be bound; stated differently, if a particle reaches infinity it is because it 
has Ut < —1. Furthermore, this condition is exact only in an axisymmetric and stationary 
spacetime, and our spacetimes attain these properties only in the final stages of the evolution. 
Nevertheless, this is a useful condition for a first estimate of the amount of matter ejected 
and a in-depth discussion on the assumptions implicit in this criterion and on how it applies 
if one accounts for external forces are presented in ||72|. (Note that the alternative criterion 
for bound flows, namely hut > — 1, would yield similar results since in the relevant regions 
/i- 1.) 

Figure [13] shows the evolution of ut embedded in a spacetime diagram much like the 
ones presented before for the rest-mass density, the specific angular momentum and the 
angular velocity. For all models under consideration, the criterion ut > —1 is well fulfilled, 
namely all the matter in the tori is bounded, except for model M3 . 4q0 . 80 which clearly 
shows in the early stages of its evolution, that a certain amount of unbound matter is ejected 
before reaching the regime of QSA. Only for the outermost, very low-density regions of the 
tori (which are not shown in the spacetime diagrams of figure [T3j values of wt < —1 are 
encountered in the other models and are probably the manifestation of an outflowing wind 
caused by the very large temperatures of those regions. As a final remark we note that 
although the total amount of matter ejected in this way is rather small and only of the order 
of - 10-"^ M0,it can nevertheless act as the site for the production of the neutron-rich heavy 
elements that are formed by rapid neutron capture [i.e. the r-process) (see [73] and references 
therein). Performing such calculations and thus determining to what extent binary NS mergers 
contribute to the whole observed r-process material in the Galaxy requires a fully developed 
reaction network and is outside the scope of this study, but will be the focus of our future 
research. 

4.8. A phenomenological expression for the mass in the torus 

As mentioned above, determining the amount of rest-mass in the torus may be one of the most 
important aspects of this research for the impact it has on the modelling of the emission in 
SGRBs. Table 12] also reports the mass of the torus and since the latter slightly decreases in 
time, we have arbitrarily chosen the time of t 17 ms as a reference (it is about the latest 
time for which we have data from all the simulations). Because of the importance of the 
information and because of the scarcity of the numerical data available, it would be valuable 
to derive a phenomenological expression for the mass in the torus which can be constructed 
on basic expectations and that can be constrained by using the numerical data. 

Following this spirit, we first search for a phenomenological expression for the torus 
mass which will depend only on the mass ratio and on the total mass of the binary, i.e. A/toi = 
Mio^{q, Aftot)- Next, we exclude the trivial case in which the total mass is larger than the 
maximum mass of the binary system A/,„ax (based on the maximum allowed mass for isolated 
stars with the given EOS); in practice we impose that Mtor('Z, Aftot > Afmax) = for any 
value of q. Finally we impose the expectation that the mass of the torus should depend, at least 
at lowest order, on the mass ratio (this was already noted by |[T]) and yield the torus with the 
smallest possible mass for an equal-mass binary. Collecting all these constraints, our ansatz 
is 

Mtor(g, Mtot) = Ci(l - g)(Af„,ax " Mtot) + C2{M,r,^^ - A^tot) 

= [c3(l + g)A4 - A/tot] [ci(l -q)+ C2] , (6) 

where in the second expression we have written the maximum mass of the binary in terms of 
the maximum mass A/* of an isolated nonrotating star, i.e. A-Zmax = £3(1 + q)M^. 
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Note that as introduced in expression (|6]l, the coefficients C2 and C3 have a direct physical 
interpretation: C2 is proportional to the mass of the torus for equal-mass binaries, while C3 
parametrises the excess of maximum mass that can be supported in the binary because of 
the stabilizing effect produced by the nonzero spin of the stars and of the tidal potential 
(i.e. C3 is expected to be slightly larger than 1). The three coefficients, ci,C2, and C3 can 
then be computed by comparing expression (|6]l with the numerical data reported in table |2] 
as well as the one computed in ||3] for equal-mass binaries. The fitting procedure then yields 
ci = 1.115 ± 1.090, C2 = 0.039 ± 0.023, C3 = 1.139 ±0.149, with a reduced ^ 2 x 10"^ 




Figure 14. Different symbols show the torus mass Mtor measured either in the simulations 
reported here (red crosses) or in those reported in (3] (green squares). Also shown in 
the parameter space (g, A/tot) considered here is the phenomenological modelling Mtor 
suggested by expression (6). Note that to highlight the functional behaviour of the 
phenomenological expression, the x— and y— axes are shown as decreasing when moving 
to the right and to the left, respectively. 

Figure[T4] shows the torus mass Mtor either as measured in the simulations reported here 
(red crosses) or in those presented in ||3] (green squares) and against the phenomenological 
modelling Mtor suggested by expression (|6]l in the region where A/tot < A/max- Note that to 
highlight the functional behaviour of the phenomenological expression, the x— and y— axes 
are shown as decreasing when moving to the right and to the left, respectively. Overall, the 
figure shows rather generically that: 1) The mass of the torus increases with the asymmetry in 
the mass ratio; 2) Such an increase is not monotonic and for sufficiently small mass ratios the 
tidal disruption leads to tori that have a smaller mass for binaries with the same total mass; 
3) Tori with masses < 0.21 Mq have been measured and even more massive ones, i.e. with 
masses up to ^ 0.35 Mq, are possible for mass ratios q ^ 0.75 — 0.85. 

We note that somewhat similar considerations about the mass of the torus were made also 
in HI, where a different phenomenological expression for the mass of the torus was proposed. 
When applied to the data computed here, the expression suggested in U does not reproduce 
well the data and yields rather large errors. There are a number of reasons that could justify 
these differences and that are related to the different initial data chosen (ref. lH] has only 
two initial total masses which are smaller than those considered here), to the different EOSs 
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employed (ref. Ul uses cold but realistic EOSs in contrast to the ideal-fluid chosen here), and 
to the different numerical techniques adopted (ref. ||T| uses a uniform grid with rather coarse 
resolution in place of the mesh-refined grid employed here). All these differences make the 
comparison between the two calculations rather difficult, although they also motivate a closer 
comparison using at least the same initial data and the same EOSs, and which will be the 
subject of our future work. However, common conclusions of both calculations are that: The 
mass of the torus can be as large as ^ 0.1 Mq and larger; It increases with the mass asymmetry 
in the binary; It is larger for systems with smaller total mass. We believe these features are 
robust and will be also present when different initial data and EOSs are considered. 

A final note of caution must be mentioned: Although figure [T4l indicates a very good 
match between the data and expression (|6]l, it also shows that the latter is inaccurate for 
g ~ 1, where the tori masses are much smaller and the prediction leads to small but negative 
values; luckily, the regime where (|6]l is less accurate is also the least interesting one from an 
astrophysical point of view. Most importantly, however, it is clear that the attempt to produce 
a phenomenological description for the mass of the torus after having investigated only a 
small portion of the space of the parameters (especially with respect to the total mass of the 
binary) and after using as support only 8 simulations is a very demanding task and potentially 
a flawed one. However, because we believe that expression (|6]l is a reasonable description of 
the expected results, we foresee that it will reveal its robustness as additional simulations are 
performed and the coefficients will be further improved. This will indeed be the subject of 
our future work. 

5. Gravitational- Wave Emission 

Figure [TS] shows the waveforms in the two polarizations of the gravitational-wave amplitude 
ih+)22 (upper panels) and (/ix)22 (lower panels) for all the models considered and as 
computed from the gauge-invariant perturbations of a Schwarzschild spacetime. As predicted 
by the post-Newtonian approximation |,74J , the inspiral phase is characterized by harmonic 
oscillations at roughly twice the orbital frequency but that show an increase both in amplitude 
and frequency as the merger approaches. We note that the initial part of the inspiral of the 
binary M3 . 4q0 . 80 shows a comparatively larger contamination from the initial spurious 
burst of radiation. This is simply due to the fact that such a binary has been constructed with 
a comparatively larger initial violation of the constraints (i.e. the violation of the L2 norm of 
the Hamiltonian constraint is ~ 3 x 10^® and about 50% larger than the violation measured in 
other binaries). We believe that this larger initial error is also the one responsible for a longer 
time spent by this binary before the merger. 

As already discussed before, because of the very high total mass of the systems no 
transient HMNS forms, whose dynamics would have been dramatically imprinted on the 
waveforms (cf. the detailed comparison of the HMNS dynamics for different EOSs presented 
in ||3])- As a result, the post-merger waveform is essentially the one corresponding to the 
collapse of the HMNS to a BH. Indeed, as noted above, in the I0W-5 cases M3 . 5q0 . 75 and 
M3 . 4q0 . 7 0, for which a common apparent horizon is found almost simultaneously with the 
merger, the part of the waveform produced by the newly formed BH starts essentially together 
with the end of the one coming from the inspiral. 

The ringdown part of the waveform starts increasingly early for binaries with smaller 
mass ratios and its signature in the waveform is also less evident. More specifically, while 
the ringdown of the BH created after the merger can be clearly identified in the waveform 
of the equal-mass model M3 . 6ql . 00, it becomes much less clear as one scrolls down in 
the different panels of figure [TS] and it seems almost absent in model M3 . 4q0 . 7 0. Indeed 
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Figure 15. Gravitational waveforms in the two polarizations hj^ (upper panels) and h x (lower 
panels) as computed from the lowest £ = m = 2 multipole for all the binaries considered. For 
those models where it was found, the vertical dashed lines mark the time of the first detection 
of the apparent horizon. Note that as the mass ratio q decreases, the lingdown part of the 
signal starts earlier but it is also less evident because of the increasingly large accretion after 
the formation of the apparent horizon. Finally, shown as an aid to comparison, the panel of 
the binary M3.4q070 also reports with dotted lines the waveforms for the equal-mass binary 
M3 . 6ql . 00. 

it is necessary to examine M3 . 4q0 . 70 on a logarithmic scale in order to appreciate the 
presence of an exponential ringdown. We believe that this behaviour is mostly likely due to the 
copious mass accretion after the formation of the apparent horizon that becomes increasingly 
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Figure 16. Scaled power spectral densities h+{f)f^^^, for all the binaries considered when 
placed at a distance of lOOMpc. Shown also are the noise curves of the Virgo detector (dotted 
magenta), of the advanced LIGO detector (short-dashed blue) and of the planned Einstein 
Telescope (long-dashed red). 



large as the mass ratio decreases. We recall, in fact, that the mass accretion rate following 
the BH formation is highly sensitive on the mass ratio and inversely proportional to it (see 
figure|5]where this is very apparent). Under these conditions of very intense mass accretion, 
the BH is continuously "hit" by generically nonspherical flows of matter which prevent its 
natural ringdown, essentially "chocking" it. A detailed analysis on the role played by mass 
accretion on the properties of the ringdown has already been investigated in ll75l . where 
however the BH ringdown was always observed because of the intrinsically perturbative 
nature of the approach. The rather different accretion regime reached in these simulations 
suggests therefore that the dynamics observed in figure [15] reflects a nonlinear response of 
the BH that was not accessible in the work of 175)1 . Additional work is needed to clarify the 
relation between hypercritical accretion and BH ringdown and will be the subject of future 
investigations. 

A more systematic analysis of the waveforms as a function of the mass ratio is beyond 
the scope of this paper and will be considered elsewhere using more realistic or parametrised 
EOSs. Here, however, as an aid to comparison, the panel relative to the binary M3.4q070 in 
figure[T5]also reports with dotted lines the waveforms for the equal-mass binary M3 . 6ql . 
and highlights that besides the different amplitude evolution, the mass asymmetry also results 
into a different phase evolution which is likely to provide important information on the EOS. 

In addition, we show with black continuous lines in figure [16] the scaled power spectral 
densities (PSD) of i.e. h+{f)f^^'^, for all the binaries considered when placed at a 
distance of lOOMpc (see 1331 for a definition of h+{f)). Shown also here are the noise curves 
of the Virgo detector (dotted magenta), of the advanced LIGO detector l76l (short-dashed 
blue), and of the planned Einstein Telescope (ET) fTTl (long-dashed red). Since the number 
of cycles computed is very small, the peak emission is the one corresponding to the last stages 
of the inspiral, around 0.6 — 0.7 kHz for all models considered. The amplitude, however, 
depends sensitively on the mass ratio, being maximal for the high-g binaries and above the 
noise curve for Virgo in these cases. As the mass ratio is decreased, in fact, the peak values 
of the PSD decrease and the binaries at the distances considered become then undetectable by 
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an interferometer like Virgo (we recall that the binary M3 . 4q0 . 8 has an extended inspiral 
induced by the larger initial violation of the constraints; hence its PSD amplitude is spuriously 
increased in figure[T6b. New-generation detectors such as advanced LIGO will instead be able 
to reveal the inspiral signal in the frequency interval ~' 0.3 — 2.0 kHz, while essentially all of 
the late-inspiral and merger signal would be measured by the Einstein Telescope (see ref. 1781 
for an introduction to the science reach of third-generation detectors such as ET, and ref 1791 
for a detailed discussion of the impact that gravitational waves from NSs may have on such 
detectors). Table |3] summarizes this by reporting the signal-to-noise ratios (SNRs) for the 
different detectors and clearly highlights that present detectors are unlikely to detect any of 
the binaries considered here if at a distance of 100 Mpc and observed only during the final 
part of the inspiral. On the other hand, advanced detectors will be able to reveal these sources 
even at such large distances (the only ones that can provide an interesting event rate) and, in 
the case of third-generation detectors such as ET, even measure them with significant SNRs. 



Model 




SNR for Virgo 


SNR for adLIGO 


SNR for ET 


M3 . 6ql . 


.00 


0.41 


2.56 


47.47 


MS . 7q0 . 


, 94 


0.41 


2.59 


48.33 


M3 . 4q0 . 


. 91 


0.38 


2.48 


45.40 


M3 . 4q0 . 


.80 


0.46 


3.29 


55.68 


M3 . 5q0 . 


, 75 


0.36 


2.38 


42.56 


M3 . 4q0 . 


.70 


0.34 


2.29 


40.48 



Table 3. SNR for the different binaries considered as computed wlien placed at a distance of 
lOOMpc for a presently operating detector such as Virgo, as well as for detectors of second 
and third generation, such as advanced LIGO and ET. 



6. Conclusions 

Numerical-relativity simulations of non-vacuum spacetimes have now reached a sufficient 
stability and accuracy to be able to describe in a complete manner all of the stages of the 
inspiral, merger and post-merger of binary NSs. Determining the properties of the black-hole- 
torus system produced by the merger represents a key aspect in the modelling of the central 
engine of SGRBs. Of the many different properties characterizing the torus, the total rest-mass 
clearly represents the most important one, since it is the torus' binding energy which can be 
tapped to extract the large energies necessary to power the SGRB emission. However, the 
rest-mass density and angular momentum distributions in the torus also represent important 
elements which determine its secular evolution and need to be computed equally accurately 
for any satisfactory modelling of SGRB engine. 

As a first step towards modelling ab-initio the central engine of SGRBs, we have 
here presented new results from accurate and fully general-relativistic simulations of the 
coalescence of unmagnetized binary NSs with unequal masses as these are the ones expected 
to yield the largest tori. The evolution of the stars has been followed through the inspiral 
phase, the merger and prompt collapse to a BH, up until the appearance of a thick accretion 
disk, which was studied as it enters and remains in a regime of quasi-steady accretion. 
Although we have employed a simple ideal-fluid equation of state, we have performed a 
systematic study of the properties of the black-hole-torus obtaining a number of results that 
can be summarized as follows: 
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• The mass of the torus increases considerably with the mass asymmetry and equal-mass 
binaries do not produce significant tori if they have a total baryonic mass Aftot ^ 
3.7 Mq. Those produced have masses Mtor ^ 10"'^ Mq and a radial extension of 
- 30 km. 

• Tori with masses as large as ~ 0.2 AIq have been measured with binaries having 
A/tot ^ 3.4 A/0 and mass ratios q ^ 0.75 — 0.85. The tori in these cases are much 
more extended with typical sizes > 120 km. 

• The mass of the torus can be described accurately by the simple expression 

Mtoriq,Mtot) = [c3(l + g)A/* - A/tot] [ci(l - + C2], involving the maximum mass 
for the binaries and some coefficients, both of which can be constrained from the 
simulations. 

• Using the phenomenological expression we conclude that tori with masses as large as 
Mtoi ~ 0.35 Mq can be produced for binaries with total masses A/tot ^ 2.8 A/q and 
g~ 0.75 -0.85. 

• Tori from equal-mass binaries exhibit a quasi-periodic form of accretion associated with 
the radial epicycUc oscillations of the tori, while those from equal-mass binaries exhibit 
a quasi-steady form of accretion. 

• When analyzing the evolution of the angular-momentum distribution in the tori, we find 
no evidence for the onset of non-axisymmetric instabilities, that angular momentum is 
transported outwards more efficiently for smaller values of q thus yielding Keplerian 
angular-velocity distributions, and that very little of the mass of the tori is unbound. 

• Present gravitational-wave detectors are unlikely to detect any of the binaries considered 
here if at a distance of 100 Mpc and observed only during the final part of the inspiral. 

• Advanced detectors will be able to reveal these sources even at large distances and 
measure them with significant SNRs in the case of third-generation detectors such as 
ET 

Overall, these results indicate that large-scale tori with large masses and quasi-stationary 
evolutions can be produced as the result of the inspiral and merger of binary NSs with unequal 
masses. Hence, they may provide the energy reservoir needed to power short GRBs. Although 
complete and accurate, our results are also far from being realistic. Much remains to be done 
to improve them either by considering physically-motivated EOSs, or by including the effect 
of magnetic fields, or by taking into account the modifications introduced by a self-consistent 
treatment of the radiation transfer. All and each of these improvements will be the subject of 
our future research. 
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Appendix A. On our accuracy: conservation of mass and angular momentum 

In a recent work jS] we have discussed in detail the convergence properties of our numerical 
simulations and, in particular, the deterioration of the convergence rate at the merger and 
during the survival of the merged object, when strong shocks are formed and turbulence 
develops. In particular, in figure 3 of that work we have shown a very stringent measure 
of the overall conservation properties of our simulations by reporting the time evolution of the 
energy and angular momentum which are partially radiated during the simulation. In order to 
reduce the computational costs associated with the measurements made in |5J, we had limited 
ourselves to a single configuration and in particular one that, because of the rather high mass, 
formed a BH soon after the merger. In particular, we had considered an equal-mass binary 
with a total baryonic mass of Mi, — 3.56 Mq and a total ADM mass of M^^^ — 3.23 Mq as 
evolved from an initial (coordinate) separation of ^ 45 km. As a result, we were able to show 
an overall conservation of both mass and angular momentum to a precision of '-^ 1% over a 
timescale of ~ 10ms. 




Figure Al. Left panel: Evolution of the maximum rest-mass density normalized to its initial 
value. Shown with different colours are the different parts of the evolution which are then also 
magnified in the three insets (cf. the timescale to associate the insets to the different parts of 
the evolution). Right panel: The same as in the left panel but in terms of the t = m = 2 mode 
of the ft+ polarization amplitude. 

We here reconsider the same assessment of the conservation properties but for a far more 
challenging case of a binary with a small total mass and for which the HMNS survives a 
considerably larger time before collapsing to a BH. In particular, we examine the evolution 
of an equal-mass binary with a total baryonic mass of A/t, = 2.912 M0 and a total ADM 
mass of A/adm ^ 2.694 M0 as evolved from an initial (coordinate) separation of ^ 45 km; 
this same binary was indicated as 1.46-45-IF in Q and evolved there only up to 25 ms. As 
representative information about the binary inspiral and merger we report in the left panel 
of figure IaTI the evolution of the maximum rest- mass density normalized to its initial value 
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{cf. figure 15 in |[3l)- Shown with different colours are the different parts of the evolution and 
which are also magnified in the different insets {cf. the timescale to associate the insets to the 
three parts of the evolution). They refer to the immediate formation of the HMNS (red line), 
to the secular evolution of the HMNS as a contracting bar-deformed object (green line) and 
to the exponential growth when the threshold to black-hole formation has been crossed (blue 
line). The right panel of the same figure shows instead the same stages of the evolution but in 
terms of the £ — m ~ 2 amplitude of the + polarization. Note that the timescale over which 
the evolution is reported is ^ 140 ms and thus a factor ^ 5 larger than the one discussed 
in |[3l. (As for the evolutions in Q, here too we have used a rotational symmetry around the 
z-axis to reduce computational costs.) 
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Figure A2. Left panel: conservation of energy. The continuous black line is the ADM mass 
computed as an integral over the whole grid, while the long-dashed blue line is the energy 
carried from gravitational waves outside the grid and magnified by a factor of 10; the short- 
dashed red line is the sum of the two and it should be conserved. The numerical violation 
is at most 1.5% (cf. dot-dashed line). Right panel: the same as in the left panel but for the 
conservation of the angular momentum. Also in this case the violation is at most 1%. 

In analogy with figure 3 of ||5l, we show in the left panel of figure lA2l the evolution of 
the total mass as normalized to the initial value (cf. left panel of figure 3 of Q). Indicated 
with different lines are the volume-integrated values of the ADM mass (solid black line), of 
the energy lost to gravitational waves magnified of a factor 10 (long-dashed blue line), and of 
their sum (short-dashed red line). The last quantity should be strictly constant and this is the 
case to a precision of ^ 0.5% during the inspiral, but with a secular decrease that brings the 
total error to be ~ 1 . 5% at the end of the simulation (as an aid to comparison the value at 1 . 1 5 
is shown with a dot-dashed line). Similar considerations apply also to the conservation of the 
angular momentum as shown in the right panel of figure lA2l (c/: right panel of figure 3 of Q) 
which uses the same conventions as the left panel (here Jvoi is computed with the integral 
(15) in Q). In this case the radiative losses are much larger (almost 15% of the available 
angular momentum is lost to gravitational waves) but the overall conservation is still accurate 
to ~ 1%. Once again it is worth nothing that the timescale over which we can show accurate 
conservation of mass and angular momentum is a factor ^ 14 larger than the one discussed 
in Is] and provides us with great confidence over the numerical accuracy of our results. Of 
course this does not provide us with any measure of whether such results are indeed realistic. 
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